Testing batch effects
#remove samples with 0 reads
ps.batch = prune_samples(sample_sums(ps.b)>=1, ps.b)
readsumsdf_samples<-data.frame(nreads = sample_sums(ps.batch),
samples = rownames(as.data.frame(sample_sums(ps.batch))),
batch = sample_data(ps.batch)$Batch.no.)
readsumsdf_samples_f = readsumsdf_samples[- grep("Control", readsumsdf_samples$batch),]
my_comparisons=list(c("1","2"),c("1","3"),c("1","4"),c("1","5"),
c("2","3"), c("2","4"), c("2","5"),
c("3","4"), c("3","5"),
c("4","5"))
ggplot(readsumsdf_samples_f, aes(x=batch, y=nreads,
color=batch, fill=batch)) +
geom_violin(alpha=0.3) +
geom_jitter(shape=16, position = "jitter") +
stat_compare_means(comparisons = my_comparisons,
method = "t.test",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Batch Number") +
guides(fill=FALSE) +
xlab("Extraction batch") +
ylab("Number of reads obtained") +
theme(axis.text.x = element_text(vjust = 1,
hjust = 1,
size = 12),
strip.text = element_text(size = 12, margin = margin()))

2 - Evaluating dataset contamination with decontam
#remove samples with 0 reads
ps.b = prune_samples(sample_sums(ps.b)>=1, ps.b)
sample_data(ps.b)$is.neg <- sample_data(ps.b)$Sample_type == "CONTROL"
contamdf0.5.prev <- isContaminant(ps.b, method="prevalence", threshold = 0.1,
neg="is.neg")
ps.final.noncontam <- prune_taxa(!contamdf0.5.prev$contaminant, ps.b)
ps.final.noncontam.noconts = subset_samples(ps.final.noncontam, Sample_type != "CONTROL")
ps.final.noncontam.noconts.nodoubs = subset_samples(ps.final.noncontam.noconts,
Sample_ID != "D02-DNA-2" & Sample_ID != "B11-DNA-7" & Sample_ID != "C10-cDNA-34" & Sample_ID != "A08-cDNA-6")
ps.final.noncontam.noconts.nodoubs = prune_samples(sample_sums(ps.final.noncontam.noconts.nodoubs)>=1000,
ps.final.noncontam.noconts.nodoubs)
ps.final <- prune_taxa(taxa_sums(ps.final.noncontam.noconts.nodoubs) > 0, ps.final.noncontam.noconts.nodoubs)
ps.final = prune_samples(sample_sums(ps.final)>=1, ps.final)
ps.final.r <- transform_sample_counts(ps.final, function(x) x / sum(x))
ps.final.r.f = filter_taxa(ps.final.r, function(x) sum(x) > .001, TRUE)
ps.final.r
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 35534 taxa and 65 samples ]
## sample_data() Sample Data: [ 65 samples by 47 sample variables ]
## tax_table() Taxonomy Table: [ 35534 taxa by 6 taxonomic ranks ]
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 2886 taxa and 65 samples ]
## sample_data() Sample Data: [ 65 samples by 47 sample variables ]
## tax_table() Taxonomy Table: [ 2886 taxa by 6 taxonomic ranks ]
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 35534 taxa and 65 samples ]
## sample_data() Sample Data: [ 65 samples by 47 sample variables ]
## tax_table() Taxonomy Table: [ 35534 taxa by 6 taxonomic ranks ]
sum(rowSums(otu_table(ps.final)))
## [1] 6750655
gl_ice_spls = as.character(unique(sample_data(ps.final)[str_detect(sample_data(ps.final)$Sample_Name, "Gl ice"),]$Sample_Name))
Setting order for all plots
sample_data(ps.final.r)$Month = factor(sample_data(ps.final.r)$Month,
levels = c("April","June","Early July", "Late July"))
sample_data(ps.final.r)$Fox_Aspect = factor(sample_data(ps.final.r)$Fox_Aspect,
levels = c("NW_AWS","N_NE","SWS1SE"))
sample_data(ps.final.r)$Snowpack_Layer = factor(sample_data(ps.final.r)$Snowpack_Layer,
levels = c("TOP","MID","BASAL","SUP ICE","GL ICE"))
sample_data(ps.final)$Month = factor(sample_data(ps.final)$Month,
levels = c("April","June","Early July", "Late July"))
sample_data(ps.final)$Fox_Aspect = factor(sample_data(ps.final)$Fox_Aspect,
levels = c("NW_AWS","N_NE","SWS1SE"))
sample_data(ps.final)$Snowpack_Layer = factor(sample_data(ps.final)$Snowpack_Layer,
levels = c("TOP","MID","BASAL","SUP ICE","GL ICE"))
sample_data(ps.final.r.f)$Month = factor(sample_data(ps.final.r.f)$Month,
levels = c("April","June","Early July", "Late July"))
sample_data(ps.final.r.f)$Fox_Aspect = factor(sample_data(ps.final.r.f)$Fox_Aspect,
levels = c("NW_AWS","N_NE","SWS1SE"))
sample_data(ps.final.r.f)$Snowpack_Layer = factor(sample_data(ps.final.r.f)$Snowpack_Layer,
levels = c("TOP","MID","BASAL","SUP ICE","GL ICE"))
3a - Phyla- and Proteobacterial class-level stacked barchart by DNA_or_cDNA~Fox_Aspect, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed
ps.final.r.glom <- tax_glom(subset_samples(ps.final.r,
Snowpack_Layer != "GL ICE" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
!(Sample_Name %in% gl_ice_spls)), taxrank = 'Phylum', NArm = FALSE)
ps.final.r.glom.psdf <- data.table(psmelt(ps.final.r.glom))
ps.final.r.glom.psdf$Phylum <- as.character(ps.final.r.glom.psdf$Phylum)
ps.final.r.glom.psdf[, mean := mean(Abundance, na.rm = TRUE),
by = "Phylum"]
ps.final.r.glom.psdf[(mean <= 0.01), Phylum := "Other"]
ps.final.r.glom.psdf$Phylum[is.na(ps.final.r.glom.psdf$Phylum) & ps.final.r.glom.psdf$Kingdom=="Bacteria"] <- "Unc. Bacteria"
ps.final.r.glom.psdf$Phylum[is.na(ps.final.r.glom.psdf$Phylum) & ps.final.r.glom.psdf$Kingdom=="Archaea"] <- "Unc. Archaea"
#Removing Proteobacteria from previous table
ps.final.r.glom.psdf.noProteo<-ps.final.r.glom.psdf[!grepl("Proteobacteria", ps.final.r.glom.psdf$Phylum),]
#New table, with Proteobacterial classes only
ps.final.r.ProtCl = subset_taxa(ps.final.r, Phylum == "Proteobacteria")
ps.final.r.ProtClglom <- tax_glom(ps.final.r.ProtCl, taxrank = 'Class', NArm = FALSE)
ps.final.r.ProtClglom.psdf <- data.table(psmelt(ps.final.r.ProtClglom))
ps.final.r.ProtClglom.psdf$Class <- as.character(ps.final.r.ProtClglom.psdf$Class)
ps.final.r.ProtClglom.psdf[, mean := mean(Abundance, na.rm = TRUE),
by = "Class"]
ps.final.r.ProtClglom.psdf[(mean <= 0.001), Class := "Other Proteobacteria"]
ps.final.r.ProtClglom.psdf.noP <- subset(ps.final.r.ProtClglom.psdf, select = -Phylum)
names(ps.final.r.glom.psdf.noProteo)[names(ps.final.r.glom.psdf.noProteo) == 'Phylum'] <- 'Taxa'
names(ps.final.r.ProtClglom.psdf.noP)[names(ps.final.r.ProtClglom.psdf.noP) == 'Class'] <- 'Taxa'
# ps.final.r.glom.psdf.noProteo = subset(ps.final.r.glom.psdf.noProteo, select=-c(Class))
ps.final.r.ProteoClAllPhy <- rbind(ps.final.r.glom.psdf.noProteo, ps.final.r.ProtClglom.psdf.noP)
tol18rainbow=c("#771155", "#AA4488", "#CC99BB", "#114477", "#4477AA", "#77AADD", "#117777", "#44AAAA", "#77CCCC", "#777711", "#AAAA44", "#DDDD77", "#774411", "#AA7744", "#DDAA77", "#771122", "#AA4455")
ps.final.r.ProteoClAllPhy$Taxa = factor(ps.final.r.ProteoClAllPhy$Taxa,
levels = c("Alphaproteobacteria", "Deltaproteobacteria",
"Gammaproteobacteria",
"Other Proteobacteria",
"Actinobacteria",
"Acidobacteria","Bacteroidetes",
"Chloroflexi","Cyanobacteria",
"Firmicutes","Gemmatimonadetes","Other"))
ggplot(ps.final.r.ProteoClAllPhy,
aes(x = Month, y = Abundance, fill = Taxa)) +
facet_grid(DNA_or_cDNA~Fox_Aspect,
scales = "free",
space = "free") +
geom_bar(stat = "identity",
position = "fill") +
scale_fill_manual(values = tol18rainbow,
na.value = "gray40") +
scale_y_continuous(expand = c(0,0),
labels = scales::percent,
breaks = c(0.25,0.5,0.75,1)) +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(angle = 35, hjust = 1),
strip.text.x = element_text(face = "bold",
size = 12),
strip.background = element_rect(size=1),
legend.text = element_text(size= 12)) +
guides(fill = guide_legend(keywidth = 1, keyheight = 1)) +
ylab("Relative Abundance\n")

3c - Phyla- and Proteobacterial class-level stacked barchart by DNA_or_cDNA~Month
tol17rainbow=c("#771155", "#AA4488", "#CC99BB", "#114477", "#4477AA", "#77AADD", "#117777", "#44AAAA", "#77CCCC", "#777711", "#AAAA44", "#DDDD77", "#774411", "#AA7744", "#DDAA77", "#AA4455")
ggplot(ps.final.r.ProteoClAllPhy,
aes(x = Sample_Name, y = Abundance, fill = Taxa)) +
facet_grid(DNA_or_cDNA~Month,
scales = "free",
space = "free") +
geom_bar(stat = "identity",
position = "fill",
color="black") +
scale_fill_manual(values = tol17rainbow,
na.value = "gray40") +
scale_y_continuous(expand = c(0,0),
labels = scales::percent,
breaks = c(0.25,0.5,0.75,1)) +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(angle = 35, hjust = 1),
strip.text.x = element_text(face = "bold",
size = 16),
strip.background = element_rect(size=1),
legend.text = element_text(size= 14)) +
guides(fill = guide_legend(keywidth = 1, keyheight = 1,
nrow = )) +
ylab("Relative Abundance\n")

5b - Genus-level stacked barchart by DNA_or_cDNA~Snowpack_Layer (filter)
Oh good. But this need improving. Getting the unclassified component of the families for the top main taxa in it, but without a guarantee for them popping into it. Also, basing it on the >0.1% dataset as per Arwyn’s advice.
ps.final.r.glom.genus.psdf_2 <- data.table(speedyseq::psmelt(ps.final.r.glom.genus))
ps.final.r.glom.genus.psdf_2$Genus <- as.character(ps.final.r.glom.genus.psdf_2$Genus)
ps.final.r.glom.genus.psdf_2[, mean := mean(Abundance, na.rm = TRUE),
by = "Genus"]
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Sphingomonadaceae"] <- "Unc. Sphingomonadaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Burkholderiaceae"] <- "Unc. Burkholderiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Labraceae"] <- "Unc. Labraceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Acetobacteraceae"] <- "Unc. Acetobacteraceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Microbacteriaceae"] <- "Unc. Microbacteriaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Deinococcaceae"] <- "Unc. Deinococcaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Gemmatimonaceae"] <- "Unc. Gemmatimonaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Intrasporangiaceae"] <- "Unc. Intrasporangiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Nocardiaceae"] <- "Unc. Nocardiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Ktedonobacteraceae"] <- "Unc. Ktedonobacteraceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Cellulomonadaceae"] <- "Unc. Cellulomonadaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Beijerinckiaceae"] <- "Unc. Beijerinckiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & is.na(ps.final.r.glom.genus.psdf_2$Family)] <- "Unclassified"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & is.na(ps.final.r.glom.genus.psdf_2$Family) & ps.final.r.glom.genus.psdf_2$Order=="Frankiales"] <- "Unc. Frankiales"
ps.final.r.glom.genus.psdf_2[(mean <= 0.01), Genus := "Other (<1%)"]
ps.final.r.glom.genus.psdf_2$Genus = factor(ps.final.r.glom.genus.psdf_2$Genus,
levels = c("Acidiphilium","Arthrobacter",
"Bacillus",
"Bosea",
"Burkholderia-Caballeronia-Paraburkholderia",
"Cryobacterium","Cupriavidus",
"Delftia","Hymenobacter",
"Labrys","Leptothrix",
"Massilia","Methylobacterium",
"Noviherbaspirillum","Oryzihumus",
"Polymorphobacter",
"Salinibacterium",
"Sphingomonas","Other (<1%)"))
tol17rainbow=c("#771155", "#AA4488", "#CC99BB", "#114477",
"#4477AA", "#77AADD", "#117777", "#44AAAA",
"#77CCCC", "#777711", "#AAAA44", "#DDDD77",
"#774411", "#AA7744", "#DDAA77", "#771122",
"#AA4455", "#DD7788", "azure3")
ggplot(ps.final.r.glom.genus.psdf_2,
aes(x = Month, y = Abundance, fill = Genus)) +
facet_grid(DNA_or_cDNA~Snowpack_Layer,
scales = "free",
space = "free") +
geom_bar(stat = "identity",
position = "fill",
color="black") +
scale_fill_manual(values = tol17rainbow,
na.value = "gray40") +
scale_y_continuous(expand = c(0,0),
labels = scales::percent,
breaks = c(0.25,0.5,0.75,1)) +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(angle = 35, hjust = 1),
strip.text.x = element_text(face = "bold",
size = 16),
strip.background = element_rect(size=1),
legend.text = element_text(size= 14)) +
guides(fill = guide_legend(keywidth = 1, keyheight = 1,
nrow = )) +
ylab("Relative Abundance\n")

5c - Genus-level stacked barchart by DNA_or_cDNA~Fox_Aspect (filter)
ggplot(ps.final.r.glom.genus.psdf_2,
aes(x = Month, y = Abundance, fill = Genus)) +
facet_grid(DNA_or_cDNA~Fox_Aspect,
scales = "free",
space = "free") +
geom_bar(stat = "identity",
position = "fill",
color="black") +
scale_fill_manual(values = tol17rainbow,
na.value = "gray40") +
scale_y_continuous(expand = c(0,0),
labels = scales::percent,
breaks = c(0.25,0.5,0.75,1)) +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(angle = 35, hjust = 1),
strip.text.x = element_text(face = "bold",
size = 16),
strip.background = element_rect(size=1),
legend.text = element_text(size= 14)) +
guides(fill = guide_legend(keywidth = 1, keyheight = 1,
nrow = )) +
ylab("Relative Abundance\n")

Lest us not forget that size of Other (<1%) corresponds to the different samples included in that category.
I don’t like this, but plotting all samples now, without grouping.
# ps.final.r.glom.genus.psdf_2$Sample_Name <- as.character(ps.final.r.glom.genus.psdf_2$Sample_Name)
# ps.final.r.glom.genus.psdf_2$Sample_Name <- factor(ps.final.r.glom.genus.psdf_2$Sample_Name,
# levels=unique(ps.final.r.glom.genus.psdf_2$Sample_Name))
ggplot(ps.final.r.glom.genus.psdf_2,
aes(x = Sample_Name, y = Abundance, fill = Genus)) +
facet_grid(DNA_or_cDNA~Month,
scales = "free",
space = "free") +
geom_bar(stat = "identity",
position = "fill",
color="black") +
scale_fill_manual(values = tol17rainbow,
na.value = "gray40") +
scale_y_continuous(expand = c(0,0),
labels = scales::percent,
breaks = c(0.25,0.5,0.75,1)) +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(angle = 35, hjust = 1),
strip.text.x = element_text(face = "bold",
size = 16),
strip.background = element_rect(size=1),
legend.text = element_text(size= 14)) +
guides(fill = guide_legend(keywidth = 1, keyheight = 1,
nrow = )) +
ylab("Relative Abundance\n")

6 - Unconstrained ordination
nmds_shapes = c(15,16,17,18,6)
nmdscolors_by_month = c("April" = "#A6048B",
"June" = "#0704A6",
"Early July" = "#727A0B",
"Late July" = "#6C0504")
nmdscolors_by_snpl = c("Top" = "#151AC6",
"Mid" = "#C6156D",
"Basal" = "#890909",
"Sup. Ice" = "#091D89",
"Glacial Ice" = "#580303")
ps.final.r.dna = subset_samples(ps.final.r, DNA_or_cDNA != "cDNA")
ps.final.r.cdna = subset_samples(ps.final.r, DNA_or_cDNA != "DNA")
ps.final.r.dna.nmds <- ordinate(ps.final.r.dna,
"NMDS",
distance="jsd",
maxit = 1e3)
## Run 0 stress 0.2107623
## Run 1 stress 0.214947
## Run 2 stress 0.2149514
## Run 3 stress 0.2243087
## Run 4 stress 0.2176892
## Run 5 stress 0.2118836
## Run 6 stress 0.2034066
## ... New best solution
## ... Procrustes: rmse 0.1219555 max resid 0.4302237
## Run 7 stress 0.2150477
## Run 8 stress 0.2102508
## Run 9 stress 0.2091659
## Run 10 stress 0.2107622
## Run 11 stress 0.2171913
## Run 12 stress 0.2164589
## Run 13 stress 0.2098655
## Run 14 stress 0.2034062
## ... New best solution
## ... Procrustes: rmse 0.0002446002 max resid 0.001304558
## ... Similar to previous best
## Run 15 stress 0.2126387
## Run 16 stress 0.228024
## Run 17 stress 0.2254012
## Run 18 stress 0.212504
## Run 19 stress 0.2184797
## Run 20 stress 0.2179839
## *** Solution reached
ps.final.r.cdna.nmds <- ordinate(ps.final.r.cdna,
"NMDS",
distance="jsd",
maxit = 1e3)
## Run 0 stress 0.1654426
## Run 1 stress 0.1628229
## ... New best solution
## ... Procrustes: rmse 0.04769025 max resid 0.162577
## Run 2 stress 0.1641699
## Run 3 stress 0.1653273
## Run 4 stress 0.1724367
## Run 5 stress 0.1659068
## Run 6 stress 0.1640257
## Run 7 stress 0.1649862
## Run 8 stress 0.1604191
## ... New best solution
## ... Procrustes: rmse 0.0849541 max resid 0.3370751
## Run 9 stress 0.1713025
## Run 10 stress 0.166409
## Run 11 stress 0.1681634
## Run 12 stress 0.1650155
## Run 13 stress 0.1754389
## Run 14 stress 0.1598542
## ... New best solution
## ... Procrustes: rmse 0.06484638 max resid 0.2497052
## Run 15 stress 0.1740954
## Run 16 stress 0.1622533
## Run 17 stress 0.1643391
## Run 18 stress 0.1653448
## Run 19 stress 0.1609927
## Run 20 stress 0.1654456
## *** No convergence -- monoMDS stopping criteria:
## 19: stress ratio > sratmax
## 1: scale factor of the gradient < sfgrmin
ps.final.r.dna.nmds.plot <- plot_ordination(ps.final.r.dna, ps.final.r.dna.nmds) +
geom_point(aes(colour=Month, shape=Snowpack_Layer),
size=5, alpha = 0.5, stroke=1.5) +
scale_shape_manual(values=nmds_shapes) +
scale_colour_manual(values=nmdscolors_by_month) +
scale_fill_manual(values=nmdscolors_by_month) +
stat_ellipse(aes(fill=Month),
geom = "polygon",
type="t",
alpha=0.15,
linetype = 2,
show.legend = FALSE) +
ggtitle("DNA nMDS") +
guides(colour = guide_legend(title="Month"))
ps.final.r.cdna.nmds.plot <- plot_ordination(ps.final.r.cdna, ps.final.r.cdna.nmds) +
geom_point(aes(colour=Month, shape=Snowpack_Layer),
size=5, alpha = 0.7, stroke=1.5) +
scale_shape_manual(values=nmds_shapes) +
scale_colour_manual(values=nmdscolors_by_month) +
scale_fill_manual(values=nmdscolors_by_month) +
stat_ellipse(aes(fill=Month),
geom = "polygon",
type="t",
alpha=0.2,
linetype = 2,
show.legend = FALSE)+
ggtitle("cDNA nMDS") +
guides(colour = guide_legend(title="Month"))
plot_grid(ps.final.r.dna.nmds.plot + theme_bw(),
ps.final.r.cdna.nmds.plot + theme_bw())

7 - Alpha-diversity metrics.
month_comparisons = list( c("April", "June"), c("April", "Early July"),
c("April", "Late July"),
c("June", "Early July"), c("June", "Late July"),
c("Early July", "Late July"))
site_comparisons = list( c("NW_AWS", "N_NE"),
c("NW_AWS", "SWS1SE"),
c("N_NE", "SWS1SE"))
slayer_comparisons = list( c("TOP", "MID"), c("TOP", "BASAL"), c("TOP", "SUP ICE"),
c("TOP", "GL ICE"),
c("MID", "BASAL"), c("MID", "SUP ICE"),
c("MID", "GL ICE"),
c("BASAL", "SUP ICE"), c("BASAL", "GL ICE"),
c("SUP ICE", "GL ICE"))
All measures for Sample_ID, no stats, facet by DNA/cDNA
plot_richness(ps.final,
x="Sample_ID", color="Sample_ID",
measures=c("Observed","Shannon","Simpson")) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Sample_ID") +
guides(color=guide_legend(ncol=2)) +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))

All measures by Sample_ID
mini_metadata = metadata %>%
select(c(Sample_ID, Sample_Name, Month, DNA_or_cDNA, Fox_Aspect, Snowpack_Layer))
ps.final.measures = estimate_richness(ps.final, measures = c("Observed","Shannon","Simpson")) %>%
rownames_to_column(var = "Sample_ID") %>%
mutate(Sample_ID = gsub("\\.", "-", Sample_ID))
full_join(ps.final.measures, mini_metadata,
by="Sample_ID") %>%
filter(!is.na(Observed)) %>%
kable(digits = 4) %>%
kable_styling(full_width = T, bootstrap_options = c("striped", "hover")) %>%
row_spec(0, bold=TRUE)
|
Sample_ID
|
Observed
|
Shannon
|
Simpson
|
Sample_Name
|
Month
|
DNA_or_cDNA
|
Fox_Aspect
|
Snowpack_Layer
|
|
A03-cDNA-1
|
286
|
0.7114
|
0.2103
|
April_N_NE_Basal_R
|
April
|
cDNA
|
N_NE
|
BASAL
|
|
A03-DNA-1
|
138
|
3.2774
|
0.9407
|
June_N_NE_Top20_D
|
June
|
DNA
|
N_NE
|
TOP
|
|
A04-cDNA-2
|
163
|
2.0361
|
0.7852
|
April_NW_AWS_Top_R
|
April
|
cDNA
|
NW_AWS
|
TOP
|
|
A04-DNA-2
|
75
|
1.6409
|
0.6083
|
April_ NW_AWS_ Basal_D
|
April
|
DNA
|
NW_AWS
|
BASAL
|
|
A05-cDNA-3
|
515
|
2.2067
|
0.7048
|
June_NW_AWS_Basal_R
|
June
|
cDNA
|
NW_AWS
|
BASAL
|
|
A05-DNA-3
|
1390
|
5.1230
|
0.9813
|
30 July_ N_NE_ Gl ice_D
|
Late July
|
DNA
|
N_NE
|
GL ICE
|
|
A06-cDNA-4
|
125
|
1.3510
|
0.6058
|
June_N_NE_Basal_R
|
June
|
cDNA
|
N_NE
|
BASAL
|
|
A06-DNA-4
|
210
|
3.9337
|
0.9650
|
June_ SWS1SE_ Top20_D
|
June
|
DNA
|
SWS1SE
|
TOP
|
|
A07-cDNA-5
|
1635
|
4.6591
|
0.9488
|
9 July_NW_AWS_Mid_R
|
Early July
|
cDNA
|
NW_AWS
|
MID
|
|
A07-DNA-5
|
609
|
3.3602
|
0.8726
|
9 July_ NW_AWS_ Sup ice_D
|
Early July
|
DNA
|
NW_AWS
|
SUP ICE
|
|
A08-DNA-6
|
563
|
4.0141
|
0.8786
|
April_ NW_AWS_ Mid_D
|
April
|
DNA
|
NW_AWS
|
MID
|
|
A09-cDNA-7
|
1875
|
5.5902
|
0.9807
|
30 July_NW_AWS_Top20_R
|
Late July
|
cDNA
|
NW_AWS
|
TOP
|
|
A10-cDNA-8
|
842
|
4.7339
|
0.9520
|
30 July_N_NE_Gl ice_R
|
Late July
|
cDNA
|
N_NE
|
GL ICE
|
|
A10-DNA-8
|
179
|
3.3916
|
0.9493
|
April_NW_AWS_Top20_D
|
April
|
DNA
|
NW_AWS
|
TOP
|
|
A11-cDNA-9
|
3
|
0.2518
|
0.1121
|
April_NW_AWS_Mid_R
|
April
|
cDNA
|
NW_AWS
|
MID
|
|
A11-DNA-9
|
1819
|
5.3121
|
0.9764
|
9 July_SWS1SE_Gl ice 1_D
|
Early July
|
DNA
|
SWS1SE
|
GL ICE
|
|
A12-cDNA-10
|
54
|
1.4929
|
0.6133
|
June_NW_AWS_Top20_R
|
June
|
cDNA
|
NW_AWS
|
TOP
|
|
A12-DNA-10
|
2719
|
5.7179
|
0.9837
|
30 July_N_NE_Top20_D
|
Late July
|
DNA
|
N_NE
|
TOP
|
|
B01-cDNA-11
|
453
|
4.0612
|
0.9622
|
30 July_NW_AWS_Sup ice_R
|
Late July
|
cDNA
|
NW_AWS
|
SUP ICE
|
|
B01-DNA-11
|
187
|
3.5277
|
0.9546
|
June_NW_AWS_Mid_D
|
June
|
DNA
|
NW_AWS
|
MID
|
|
B02-cDNA-12
|
733
|
4.5851
|
0.9683
|
30 July_N_NE_Top20_R
|
Late July
|
cDNA
|
N_NE
|
TOP
|
|
B02-DNA-12
|
310
|
2.1151
|
0.7154
|
9 July_N_NE_Mid_D
|
Early July
|
DNA
|
N_NE
|
MID
|
|
B03-cDNA-13
|
38
|
1.1241
|
0.4091
|
June_SWS1SE_Basal_R
|
June
|
cDNA
|
SWS1SE
|
BASAL
|
|
B03-DNA-13
|
457
|
4.2859
|
0.9559
|
April_SWS1SE_Mid_D
|
April
|
DNA
|
SWS1SE
|
MID
|
|
B04-cDNA-14
|
67
|
1.1867
|
0.4139
|
9 July_NW_AWS_Top20_R
|
Early July
|
cDNA
|
NW_AWS
|
TOP
|
|
B05-cDNA-15
|
1307
|
5.7881
|
0.9881
|
9 July_NW_AWS_Gl ice 1_R
|
Early July
|
cDNA
|
NW_AWS
|
GL ICE
|
|
B06-cDNA-16
|
134
|
1.4982
|
0.6595
|
April_SWS1SE_Top20_R
|
April
|
cDNA
|
SWS1SE
|
TOP
|
|
B06-DNA-2
|
156
|
0.6586
|
0.2271
|
April_SWS1SE_Basal_D
|
April
|
DNA
|
SWS1SE
|
BASAL
|
|
B07-DNA-3
|
1302
|
5.3971
|
0.9910
|
June_NW AWS_Top20_D
|
June
|
DNA
|
NW_AWS
|
TOP
|
|
B08-cDNA-18
|
21
|
0.9250
|
0.5154
|
June_NW_AWS_Mid_R
|
June
|
cDNA
|
NW_AWS
|
MID
|
|
B08-DNA-4
|
986
|
5.2137
|
0.9805
|
30 July_SWS1SE_Gl ice_D
|
Late July
|
DNA
|
SWS1SE
|
GL ICE
|
|
B09-DNA-5
|
105
|
2.9238
|
0.9150
|
June_N_NE_Mid_D
|
June
|
DNA
|
N_NE
|
MID
|
|
B10-cDNA-21
|
400
|
4.5700
|
0.9771
|
9 July_SWS1SE_Sup ice_R
|
Early July
|
cDNA
|
SWS1SE
|
SUP ICE
|
|
B10-DNA-6
|
377
|
4.1547
|
0.9598
|
June_SWS1SE_Basal_D
|
June
|
DNA
|
SWS1SE
|
BASAL
|
|
B12-cDNA-23
|
29
|
0.9625
|
0.3511
|
30 July_SWS1SE_Slush_R
|
Late July
|
cDNA
|
SWS1SE
|
TOP
|
|
B12-DNA-8
|
1998
|
5.3382
|
0.9817
|
9 July_ NW_AWS_Gl ice 1_D
|
Early July
|
DNA
|
NW_AWS
|
GL ICE
|
|
C01-DNA-9
|
165
|
3.3370
|
0.9306
|
June_SWS1SE_Mid_D
|
June
|
DNA
|
SWS1SE
|
MID
|
|
C02-cDNA-26
|
907
|
4.1753
|
0.9419
|
30 July_SWS1SE_Gl ice_R
|
Late July
|
cDNA
|
SWS1SE
|
GL ICE
|
|
C03-cDNA-27
|
537
|
3.3393
|
0.8840
|
9 July_N_NE_Sup ice_R
|
Early July
|
cDNA
|
N_NE
|
SUP ICE
|
|
C03-DNA-B2-1
|
677
|
5.1429
|
0.9887
|
June_NW_AWS_Basal_D
|
June
|
DNA
|
NW_AWS
|
BASAL
|
|
C04-cDNA-28
|
307
|
3.3747
|
0.8749
|
30 July_NW_AWS_Gl ice_R
|
Late July
|
cDNA
|
NW_AWS
|
GL ICE
|
|
C04-DNA-2
|
715
|
4.2502
|
0.9430
|
30 July_NW_AWS_Gl ice_D
|
Late July
|
DNA
|
NW_AWS
|
GL ICE
|
|
C05-DNA-3
|
1739
|
5.2347
|
0.9838
|
30 July_NW_AWS_Sup ice_D
|
Late July
|
DNA
|
NW_AWS
|
SUP ICE
|
|
C06-DNA-4
|
236
|
3.4402
|
0.9006
|
April_SWS1SE_Top20_D
|
April
|
DNA
|
SWS1SE
|
TOP
|
|
C07-DNA-5
|
2789
|
5.9749
|
0.9894
|
9 July_NW_AWS_Mid_D
|
Early July
|
DNA
|
NW_AWS
|
MID
|
|
C08-DNA-6
|
307
|
4.5001
|
0.9837
|
June_N_NE_Basal_D
|
June
|
DNA
|
N_NE
|
BASAL
|
|
C09-cDNA-33
|
36
|
1.3284
|
0.5802
|
June_N_NE_Top20_R
|
June
|
cDNA
|
N_NE
|
TOP
|
|
C09-DNA-7
|
148
|
3.8243
|
0.9690
|
9 July_N_NE_Gl ice_D
|
Early July
|
DNA
|
N_NE
|
GL ICE
|
|
C10-DNA-8
|
383
|
3.2676
|
0.8801
|
9 July_ SWS1SE_Mid_D
|
Early July
|
DNA
|
SWS1SE
|
MID
|
|
C11-cDNA-35
|
65
|
1.5747
|
0.6657
|
9 July_SWS1SE_Mid_R
|
Early July
|
cDNA
|
SWS1SE
|
MID
|
|
C12-cDNA-36
|
19
|
0.4213
|
0.1462
|
9 July_N_NE_Gl ice_R
|
Early July
|
cDNA
|
N_NE
|
GL ICE
|
|
C12-DNA-10
|
396
|
2.2692
|
0.7492
|
9 July_N_NE_Sup ice_D
|
Early July
|
DNA
|
N_NE
|
SUP ICE
|
|
D01-DNA-B4-1
|
147
|
0.6873
|
0.2422
|
30 July_SWS1SE_Slush_D
|
Late July
|
DNA
|
SWS1SE
|
TOP
|
|
D02-cDNA-38
|
128
|
2.1030
|
0.7780
|
April_SWS1SE_Mid_R
|
April
|
cDNA
|
SWS1SE
|
MID
|
|
D03-cDNA-39
|
266
|
3.4007
|
0.9035
|
April_N_NE_Mid_R
|
April
|
cDNA
|
N_NE
|
MID
|
|
D04-cDNA-40
|
79
|
1.5083
|
0.6114
|
June_SWS1SE_Top20_R
|
June
|
cDNA
|
SWS1SE
|
TOP
|
|
D05-cDNA-41
|
154
|
1.8281
|
0.7377
|
June_N_NE_Mid_R
|
June
|
cDNA
|
N_NE
|
MID
|
|
D05-DNA-5
|
90
|
0.6634
|
0.2169
|
9 July_N_NE_Top20_D
|
Early July
|
DNA
|
N_NE
|
TOP
|
|
D06-cDNA-42
|
841
|
3.6042
|
0.9024
|
9 July_NW_AWS_Sup ice_R
|
Early July
|
cDNA
|
NW_AWS
|
SUP ICE
|
|
D06-DNA-6
|
2292
|
5.5649
|
0.9810
|
9 July_SWS1SE_Sup ice_D
|
Early July
|
DNA
|
SWS1SE
|
SUP ICE
|
|
D07-cDNA-43
|
1444
|
4.5026
|
0.9489
|
9 July_SWS1SE_Gl ice 2_R
|
Early July
|
cDNA
|
SWS1SE
|
GL ICE
|
|
D07-DNA-7
|
467
|
1.8746
|
0.5308
|
9 July_NW_AWS_Top20_D
|
Early July
|
DNA
|
NW_AWS
|
TOP
|
|
D08-cDNA-44
|
120
|
2.6011
|
0.8428
|
9 July_N_NE_Mid_R
|
Early July
|
cDNA
|
N_NE
|
MID
|
|
D08-DNA-8
|
1045
|
4.1510
|
0.9245
|
30 July_NW_AWS_Top20_D
|
Late July
|
DNA
|
NW_AWS
|
TOP
|
|
D09-DNA-9
|
2986
|
5.9104
|
0.9895
|
9 July_SWS1SE_Top20_D
|
Early July
|
DNA
|
SWS1SE
|
TOP
|
All measures for Month, stats included, facet by DNA/cDNA
plot_richness(ps.final,
x="Month", color="Month",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Month") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))

All measures for Month, stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed
plot_richness(subset_samples(ps.final,
Snowpack_Layer != "GL ICE" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
!(Sample_Name %in% gl_ice_spls)),
x="Month", color="Month",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Month") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))

Table with means
estimate_richness(subset_samples(ps.final,
Snowpack_Layer != "GL ICE" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
!(Sample_Name %in% gl_ice_spls)),
measures = c("Observed","Shannon","Simpson")) %>%
rownames_to_column(var = "Sample_ID") %>%
mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>%
full_join(., mini_metadata, by="Sample_ID") %>%
select(!c(Sample_ID, Sample_Name, Fox_Aspect, Snowpack_Layer)) %>%
filter(!is.na(Observed)) %>%
group_by(Month, DNA_or_cDNA) %>%
mutate(Observed = mean(Observed),
Shannon = mean(Shannon),
Simpson = mean(Simpson)) %>%
distinct_all() %>%
kable(digits = 4) %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>%
row_spec(0, bold=TRUE)
|
Observed
|
Shannon
|
Simpson
|
Month
|
DNA_or_cDNA
|
|
169.2000
|
1.7006
|
0.5578
|
April
|
cDNA
|
|
270.7500
|
3.8497
|
0.9548
|
June
|
DNA
|
|
277.6667
|
2.9052
|
0.7533
|
April
|
DNA
|
|
158.6000
|
1.5037
|
0.5823
|
June
|
cDNA
|
|
695.6000
|
3.5495
|
0.8756
|
Early July
|
cDNA
|
|
1279.0000
|
3.7921
|
0.8385
|
Early July
|
DNA
|
|
785.6667
|
3.5380
|
0.7647
|
Late July
|
cDNA
|
|
2229.0000
|
5.4763
|
0.9837
|
Late July
|
DNA
|
All measures for Month, stats included, facet by DNA/cDNA, no (“Early July”, “Late July”)
mini_month_comparisons = list(c("April", "June"),
c("April", "June"))
plot_richness(subset_samples(ps.final, Month != "Early July" & Month != "Late July"),
x="Month", color="Month",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = mini_month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Month") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))

Table with means
estimate_richness(subset_samples(ps.final,
Month != "Early July" & Month != "Late July"),
measures = c("Observed","Shannon","Simpson")) %>%
rownames_to_column(var = "Sample_ID") %>%
mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>%
full_join(., mini_metadata, by="Sample_ID") %>%
select(!c(Sample_ID, Sample_Name, Fox_Aspect, Snowpack_Layer)) %>%
filter(!is.na(Observed)) %>%
group_by(Month, DNA_or_cDNA) %>%
mutate(Observed = mean(Observed),
Shannon = mean(Shannon),
Simpson = mean(Simpson)) %>%
distinct_all() %>%
kable(digits = 4) %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>%
row_spec(0, bold=TRUE)
|
Observed
|
Shannon
|
Simpson
|
Month
|
DNA_or_cDNA
|
|
163.3333
|
1.6669
|
0.5748
|
April
|
cDNA
|
|
385.3333
|
4.0216
|
0.9588
|
June
|
DNA
|
|
277.6667
|
2.9052
|
0.7533
|
April
|
DNA
|
|
127.7500
|
1.4706
|
0.5972
|
June
|
cDNA
|
All measures for Fox_Aspect, no stats, facet by DNA/cDNA
#now without stats, only shapes
plot_richness(ps.final,
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(aes(shape=Month), size=5, alpha=0.7,
position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Fox Aspect") +
scale_shape_manual(values=nmds_shapes) +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))

All measures for Fox_Aspect, stats included, facet by DNA/cDNA
plot_richness(ps.final,
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha=0.5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = site_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Fox Aspect") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))

Table with means
estimate_richness(ps.final,
measures = c("Observed","Shannon","Simpson")) %>%
rownames_to_column(var = "Sample_ID") %>%
mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>%
full_join(., mini_metadata, by="Sample_ID") %>%
select(!c(Sample_ID, Sample_Name, Snowpack_Layer)) %>%
filter(!is.na(Observed)) %>%
group_by(Month, DNA_or_cDNA, Fox_Aspect) %>%
mutate(Observed = mean(Observed),
Shannon = mean(Shannon),
Simpson = mean(Simpson)) %>%
distinct_all() %>%
kable(digits = 4) %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>%
row_spec(0, bold=TRUE)
|
Observed
|
Shannon
|
Simpson
|
Month
|
DNA_or_cDNA
|
Fox_Aspect
|
|
276.0000
|
2.0560
|
0.5569
|
April
|
cDNA
|
N_NE
|
|
183.3333
|
3.5671
|
0.9465
|
June
|
DNA
|
N_NE
|
|
83.0000
|
1.1439
|
0.4487
|
April
|
cDNA
|
NW_AWS
|
|
272.3333
|
3.0155
|
0.8120
|
April
|
DNA
|
NW_AWS
|
|
196.6667
|
1.5415
|
0.6112
|
June
|
cDNA
|
NW_AWS
|
|
2054.5000
|
5.4205
|
0.9825
|
Late July
|
DNA
|
N_NE
|
|
105.0000
|
1.5025
|
0.6412
|
June
|
cDNA
|
N_NE
|
|
250.6667
|
3.8084
|
0.9518
|
June
|
DNA
|
SWS1SE
|
|
962.5000
|
3.8095
|
0.8133
|
Early July
|
cDNA
|
NW_AWS
|
|
1465.7500
|
4.1370
|
0.8436
|
Early July
|
DNA
|
NW_AWS
|
|
878.3333
|
4.3420
|
0.9393
|
Late July
|
cDNA
|
NW_AWS
|
|
787.5000
|
4.6595
|
0.9602
|
Late July
|
cDNA
|
N_NE
|
|
1870.0000
|
5.0137
|
0.9568
|
Early July
|
DNA
|
SWS1SE
|
|
722.0000
|
4.6893
|
0.9781
|
June
|
DNA
|
NW_AWS
|
|
236.0000
|
2.2180
|
0.6626
|
Early July
|
DNA
|
N_NE
|
|
58.5000
|
1.3162
|
0.5102
|
June
|
cDNA
|
SWS1SE
|
|
283.0000
|
2.7949
|
0.6945
|
April
|
DNA
|
SWS1SE
|
|
131.0000
|
1.8006
|
0.7188
|
April
|
cDNA
|
SWS1SE
|
|
566.5000
|
2.9505
|
0.6113
|
Late July
|
DNA
|
SWS1SE
|
|
636.3333
|
3.5491
|
0.8639
|
Early July
|
cDNA
|
SWS1SE
|
|
468.0000
|
2.5689
|
0.6465
|
Late July
|
cDNA
|
SWS1SE
|
|
225.3333
|
2.1206
|
0.6244
|
Early July
|
cDNA
|
N_NE
|
|
1166.3333
|
4.5453
|
0.9504
|
Late July
|
DNA
|
NW_AWS
|
All measures for Fox_Aspect, stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed
plot_richness(subset_samples(ps.final,
Snowpack_Layer != "GL ICE" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
!(Sample_Name %in% gl_ice_spls)),
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha=0.5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = site_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Fox Aspect") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))

Table with means
estimate_richness(subset_samples(ps.final,
Snowpack_Layer != "GL ICE" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
!(Sample_Name %in% gl_ice_spls)),
measures = c("Observed","Shannon","Simpson")) %>%
rownames_to_column(var = "Sample_ID") %>%
mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>%
full_join(., mini_metadata, by="Sample_ID") %>%
select(!c(Sample_ID, Sample_Name, Snowpack_Layer)) %>%
filter(!is.na(Observed)) %>%
group_by(Month, DNA_or_cDNA, Fox_Aspect) %>%
mutate(Observed = mean(Observed),
Shannon = mean(Shannon),
Simpson = mean(Simpson)) %>%
distinct_all() %>%
kable(digits = 4) %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>%
row_spec(0, bold=TRUE)
|
Observed
|
Shannon
|
Simpson
|
Month
|
DNA_or_cDNA
|
Fox_Aspect
|
|
276.0000
|
2.0560
|
0.5569
|
April
|
cDNA
|
N_NE
|
|
183.3333
|
3.5671
|
0.9465
|
June
|
DNA
|
N_NE
|
|
83.0000
|
1.1439
|
0.4487
|
April
|
cDNA
|
NW_AWS
|
|
272.3333
|
3.0155
|
0.8120
|
April
|
DNA
|
NW_AWS
|
|
515.0000
|
2.2067
|
0.7048
|
June
|
cDNA
|
NW_AWS
|
|
80.5000
|
1.3397
|
0.5930
|
June
|
cDNA
|
N_NE
|
|
250.6667
|
3.8084
|
0.9518
|
June
|
DNA
|
SWS1SE
|
|
1238.0000
|
4.1317
|
0.9256
|
Early July
|
cDNA
|
NW_AWS
|
|
1288.3333
|
3.7366
|
0.7976
|
Early July
|
DNA
|
NW_AWS
|
|
1164.0000
|
4.8257
|
0.9715
|
Late July
|
cDNA
|
NW_AWS
|
|
2719.0000
|
5.7179
|
0.9837
|
Late July
|
DNA
|
N_NE
|
|
432.0000
|
4.3353
|
0.9716
|
June
|
DNA
|
NW_AWS
|
|
353.0000
|
2.1921
|
0.7323
|
Early July
|
DNA
|
N_NE
|
|
58.5000
|
1.3162
|
0.5102
|
June
|
cDNA
|
SWS1SE
|
|
283.0000
|
2.7949
|
0.6945
|
April
|
DNA
|
SWS1SE
|
|
232.5000
|
3.0723
|
0.8214
|
Early July
|
cDNA
|
SWS1SE
|
|
29.0000
|
0.9625
|
0.3511
|
Late July
|
cDNA
|
SWS1SE
|
|
537.0000
|
3.3393
|
0.8840
|
Early July
|
cDNA
|
N_NE
|
|
1739.0000
|
5.2347
|
0.9838
|
Late July
|
DNA
|
NW_AWS
|
|
1887.0000
|
4.9143
|
0.9502
|
Early July
|
DNA
|
SWS1SE
|
|
128.0000
|
2.1030
|
0.7780
|
April
|
cDNA
|
SWS1SE
|
All measures for Fox_Aspect, stats included, facet by DNA/cDNA, no (“Early July”, “Late July”)
plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
Month != "Early July" & Month != "Late July"),
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = mini_month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Fox_Aspect") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

Table with means
estimate_richness(subset_samples(ps.final,
Snowpack_Layer != "GL ICE" &
Month != "Early July" & Month != "Late July"),
measures = c("Observed","Shannon","Simpson")) %>%
rownames_to_column(var = "Sample_ID") %>%
mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>%
full_join(., mini_metadata, by="Sample_ID") %>%
select(!c(Sample_ID, Sample_Name, Snowpack_Layer)) %>%
filter(!is.na(Observed)) %>%
group_by(Month, DNA_or_cDNA, Fox_Aspect) %>%
mutate(Observed = mean(Observed),
Shannon = mean(Shannon),
Simpson = mean(Simpson)) %>%
distinct_all() %>%
kable(digits = 4) %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>%
row_spec(0, bold=TRUE)
|
Observed
|
Shannon
|
Simpson
|
Month
|
DNA_or_cDNA
|
Fox_Aspect
|
|
276.0000
|
2.0560
|
0.5569
|
April
|
cDNA
|
N_NE
|
|
183.3333
|
3.5671
|
0.9465
|
June
|
DNA
|
N_NE
|
|
83.0000
|
1.1439
|
0.4487
|
April
|
cDNA
|
NW_AWS
|
|
272.3333
|
3.0155
|
0.8120
|
April
|
DNA
|
NW_AWS
|
|
196.6667
|
1.5415
|
0.6112
|
June
|
cDNA
|
NW_AWS
|
|
105.0000
|
1.5025
|
0.6412
|
June
|
cDNA
|
N_NE
|
|
250.6667
|
3.8084
|
0.9518
|
June
|
DNA
|
SWS1SE
|
|
722.0000
|
4.6893
|
0.9781
|
June
|
DNA
|
NW_AWS
|
|
58.5000
|
1.3162
|
0.5102
|
June
|
cDNA
|
SWS1SE
|
|
283.0000
|
2.7949
|
0.6945
|
April
|
DNA
|
SWS1SE
|
|
131.0000
|
1.8006
|
0.7188
|
April
|
cDNA
|
SWS1SE
|
All measures for Fox_Aspect, stats included, facet by DNA/cDNA, July only
plot_richness(subset_samples(ps.final, Month == c("Early July", "Late July")),
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = mini_month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Fox_Aspect") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
## Warning in `==.default`(Month, c("Early July", "Late July")): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

Table with means
estimate_richness(subset_samples(ps.final,
Month == c("Early July", "Late July")),
measures = c("Observed","Shannon","Simpson")) %>%
rownames_to_column(var = "Sample_ID") %>%
mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>%
full_join(., mini_metadata, by="Sample_ID") %>%
select(!c(Sample_ID, Sample_Name, Snowpack_Layer)) %>%
filter(!is.na(Observed)) %>%
group_by(Month, DNA_or_cDNA, Fox_Aspect) %>%
mutate(Observed = mean(Observed),
Shannon = mean(Shannon),
Simpson = mean(Simpson)) %>%
distinct_all() %>%
kable(digits = 4) %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>%
row_spec(0, bold=TRUE)
## Warning in `==.default`(Month, c("Early July", "Late July")): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
|
Observed
|
Shannon
|
Simpson
|
Month
|
DNA_or_cDNA
|
Fox_Aspect
|
|
2054.5000
|
5.4205
|
0.9825
|
Late July
|
DNA
|
N_NE
|
|
847.6667
|
3.1500
|
0.7550
|
Early July
|
cDNA
|
NW_AWS
|
|
1875.0000
|
5.5902
|
0.9807
|
Late July
|
cDNA
|
NW_AWS
|
|
922.0000
|
4.5363
|
0.9630
|
Early July
|
cDNA
|
SWS1SE
|
|
907.0000
|
4.1753
|
0.9419
|
Late July
|
cDNA
|
SWS1SE
|
|
225.3333
|
2.1206
|
0.6244
|
Early July
|
cDNA
|
N_NE
|
|
880.0000
|
4.2006
|
0.9338
|
Late July
|
DNA
|
NW_AWS
|
|
2789.0000
|
5.9749
|
0.9894
|
Early July
|
DNA
|
NW_AWS
|
|
1684.5000
|
4.5890
|
0.9348
|
Early July
|
DNA
|
SWS1SE
|
All measures for Fox_Aspect, stats included, facet by DNA/cDNA, July only, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed."
plot_richness(subset_samples(ps.final, Month == c("Early July", "Late July"),
Snowpack_Layer != "GL ICE" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
!(Sample_Name %in% gl_ice_spls)),
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = mini_month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Fox_Aspect") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
## Warning in `==.default`(Month, c("Early July", "Late July")): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

Table with means
estimate_richness(subset_samples(ps.final, Month == c("Early July", "Late July"),
Snowpack_Layer != "GL ICE" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
!(Sample_Name %in% gl_ice_spls)),
measures = c("Observed","Shannon","Simpson")) %>%
rownames_to_column(var = "Sample_ID") %>%
mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>%
full_join(., mini_metadata, by="Sample_ID") %>%
select(!c(Sample_ID, Sample_Name, Snowpack_Layer)) %>%
filter(!is.na(Observed)) %>%
group_by(Month, DNA_or_cDNA, Fox_Aspect) %>%
mutate(Observed = mean(Observed),
Shannon = mean(Shannon),
Simpson = mean(Simpson)) %>%
distinct_all() %>%
kable(digits = 4) %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>%
row_spec(0, bold=TRUE)
## Warning in `==.default`(Month, c("Early July", "Late July")): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
|
Observed
|
Shannon
|
Simpson
|
Month
|
DNA_or_cDNA
|
Fox_Aspect
|
|
2054.5000
|
5.4205
|
0.9825
|
Late July
|
DNA
|
N_NE
|
|
847.6667
|
3.1500
|
0.7550
|
Early July
|
cDNA
|
NW_AWS
|
|
1875.0000
|
5.5902
|
0.9807
|
Late July
|
cDNA
|
NW_AWS
|
|
922.0000
|
4.5363
|
0.9630
|
Early July
|
cDNA
|
SWS1SE
|
|
907.0000
|
4.1753
|
0.9419
|
Late July
|
cDNA
|
SWS1SE
|
|
225.3333
|
2.1206
|
0.6244
|
Early July
|
cDNA
|
N_NE
|
|
880.0000
|
4.2006
|
0.9338
|
Late July
|
DNA
|
NW_AWS
|
|
2789.0000
|
5.9749
|
0.9894
|
Early July
|
DNA
|
NW_AWS
|
|
1684.5000
|
4.5890
|
0.9348
|
Early July
|
DNA
|
SWS1SE
|
All measures for Fox_Aspect, stats included, facet by DNA/cDNA, “Early July” only, “GL ICE” removed, “Gl ice” samples removed
plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
Month == "Early July" & !(Sample_Name %in% gl_ice_spls)),
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = mini_month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Fox_Aspect") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

Table with means
estimate_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
Month == "Early July" & !(Sample_Name %in% gl_ice_spls)),
measures = c("Observed","Shannon","Simpson")) %>%
rownames_to_column(var = "Sample_ID") %>%
mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>%
full_join(., mini_metadata, by="Sample_ID") %>%
select(!c(Sample_ID, Sample_Name, Snowpack_Layer)) %>%
filter(!is.na(Observed)) %>%
group_by(Month, DNA_or_cDNA, Fox_Aspect) %>%
mutate(Observed = mean(Observed),
Shannon = mean(Shannon),
Simpson = mean(Simpson)) %>%
distinct_all() %>%
kable(digits = 4) %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>%
row_spec(0, bold=TRUE)
|
Observed
|
Shannon
|
Simpson
|
Month
|
DNA_or_cDNA
|
Fox_Aspect
|
|
1238.000
|
4.1317
|
0.9256
|
Early July
|
cDNA
|
NW_AWS
|
|
1288.333
|
3.7366
|
0.7976
|
Early July
|
DNA
|
NW_AWS
|
|
353.000
|
2.1921
|
0.7323
|
Early July
|
DNA
|
N_NE
|
|
232.500
|
3.0723
|
0.8214
|
Early July
|
cDNA
|
SWS1SE
|
|
537.000
|
3.3393
|
0.8840
|
Early July
|
cDNA
|
N_NE
|
|
1887.000
|
4.9143
|
0.9502
|
Early July
|
DNA
|
SWS1SE
|
All measures for Fox_Aspect, stats included, facet by DNA/cDNA, “Late July” only, “GL ICE”, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” removed, “Gl ice” samples removed
plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" & Month == "Late July" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
!(Sample_Name %in% gl_ice_spls)),
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = mini_month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Fox_Aspect") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

Table with means
estimate_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" & Month == "Late July" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
!(Sample_Name %in% gl_ice_spls)),
measures = c("Observed","Shannon","Simpson")) %>%
rownames_to_column(var = "Sample_ID") %>%
mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>%
full_join(., mini_metadata, by="Sample_ID") %>%
select(!c(Sample_ID, Sample_Name, Snowpack_Layer)) %>%
filter(!is.na(Observed)) %>%
group_by(Month, DNA_or_cDNA, Fox_Aspect) %>%
mutate(Observed = mean(Observed),
Shannon = mean(Shannon),
Simpson = mean(Simpson)) %>%
distinct_all() %>%
kable(digits = 4) %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>%
row_spec(0, bold=TRUE)
|
Observed
|
Shannon
|
Simpson
|
Month
|
DNA_or_cDNA
|
Fox_Aspect
|
|
1164
|
4.8257
|
0.9715
|
Late July
|
cDNA
|
NW_AWS
|
|
2719
|
5.7179
|
0.9837
|
Late July
|
DNA
|
N_NE
|
|
29
|
0.9625
|
0.3511
|
Late July
|
cDNA
|
SWS1SE
|
|
1739
|
5.2347
|
0.9838
|
Late July
|
DNA
|
NW_AWS
|
All measures for Snowpack_Layer, no stats, facet by DNA/cDNA
#now without stats, only shapes
plot_richness(ps.final,
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(aes(shape=Month), size=5, alpha=0.7,
position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Snowpack Layer") +
scale_shape_manual(values=nmds_shapes) +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))

All measures for Snowpack_Layer, stats included, facet by DNA/cDNA
plot_richness(ps.final,
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = slayer_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Snowpack Layer") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))

Table with means
estimate_richness(ps.final,
measures = c("Observed","Shannon","Simpson")) %>%
rownames_to_column(var = "Sample_ID") %>%
mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>%
full_join(., mini_metadata, by="Sample_ID") %>%
select(!c(Sample_ID, Sample_Name, Fox_Aspect, Month)) %>%
filter(!is.na(Observed)) %>%
group_by(DNA_or_cDNA, Snowpack_Layer) %>%
mutate(Observed = mean(Observed),
Shannon = mean(Shannon),
Simpson = mean(Simpson)) %>%
distinct_all() %>%
kable(digits = 4) %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>%
row_spec(0, bold=TRUE)
|
Observed
|
Shannon
|
Simpson
|
DNA_or_cDNA
|
Snowpack_Layer
|
|
241.0000
|
1.3483
|
0.4825
|
cDNA
|
BASAL
|
|
865.3636
|
3.4950
|
0.7849
|
DNA
|
TOP
|
|
352.2222
|
2.2432
|
0.6626
|
cDNA
|
TOP
|
|
318.4000
|
3.2194
|
0.7535
|
DNA
|
BASAL
|
|
1176.0000
|
4.8436
|
0.9720
|
DNA
|
GL ICE
|
|
299.0000
|
2.1679
|
0.6880
|
cDNA
|
MID
|
|
1259.0000
|
4.1072
|
0.8966
|
DNA
|
SUP ICE
|
|
619.8750
|
3.6808
|
0.9025
|
DNA
|
MID
|
|
804.3333
|
3.8327
|
0.8087
|
cDNA
|
GL ICE
|
|
557.7500
|
3.8937
|
0.9314
|
cDNA
|
SUP ICE
|
All measures for Snowpack_Layer in July only (“Early July”, “Late July”), stats included, facet by DNA/cDNA
plot_richness(subset_samples(ps.final, Month = c("Early July", "Late July")),
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = slayer_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Snowpack Layer") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))

# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)
Table with means
estimate_richness(subset_samples(ps.final,
Month = c("Early July", "Late July")),
measures = c("Observed","Shannon","Simpson")) %>%
rownames_to_column(var = "Sample_ID") %>%
mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>%
full_join(., mini_metadata, by="Sample_ID") %>%
select(!c(Sample_ID, Sample_Name, Fox_Aspect, Month)) %>%
filter(!is.na(Observed)) %>%
group_by(DNA_or_cDNA, Snowpack_Layer) %>%
mutate(Observed = mean(Observed),
Shannon = mean(Shannon),
Simpson = mean(Simpson)) %>%
distinct_all() %>%
kable(digits = 4) %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>%
row_spec(0, bold=TRUE)
|
Observed
|
Shannon
|
Simpson
|
DNA_or_cDNA
|
Snowpack_Layer
|
|
241.0000
|
1.3483
|
0.4825
|
cDNA
|
BASAL
|
|
865.3636
|
3.4950
|
0.7849
|
DNA
|
TOP
|
|
352.2222
|
2.2432
|
0.6626
|
cDNA
|
TOP
|
|
318.4000
|
3.2194
|
0.7535
|
DNA
|
BASAL
|
|
1176.0000
|
4.8436
|
0.9720
|
DNA
|
GL ICE
|
|
299.0000
|
2.1679
|
0.6880
|
cDNA
|
MID
|
|
1259.0000
|
4.1072
|
0.8966
|
DNA
|
SUP ICE
|
|
619.8750
|
3.6808
|
0.9025
|
DNA
|
MID
|
|
804.3333
|
3.8327
|
0.8087
|
cDNA
|
GL ICE
|
|
557.7500
|
3.8937
|
0.9314
|
cDNA
|
SUP ICE
|
All measures for Snowpack_Layer in July only (“Early July”, “Late July”), stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed, “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” removed
plot_richness(subset_samples(ps.final, Month == c("Early July", "Late July") &
Snowpack_Layer != "GL ICE" &
!(Sample_Name %in% gl_ice_spls) &
Sample_Name != c("30 July_SWS1SE_Slush_R",
"30 July_SWS1SE_Slush_D")),
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = slayer_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Snowpack Layer") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
## Warning in `==.default`(Month, c("Early July", "Late July")): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning in `!=.default`(Sample_Name, c("30 July_SWS1SE_Slush_R", "30
## July_SWS1SE_Slush_D")): longer object length is not a multiple of shorter object
## length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)
Table with means
estimate_richness(subset_samples(ps.final, Month == c("Early July", "Late July") &
Snowpack_Layer != "GL ICE" &
!(Sample_Name %in% gl_ice_spls) &
Sample_Name != c("30 July_SWS1SE_Slush_R",
"30 July_SWS1SE_Slush_D")),
measures = c("Observed","Shannon","Simpson")) %>%
rownames_to_column(var = "Sample_ID") %>%
mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>%
full_join(., mini_metadata, by="Sample_ID") %>%
select(!c(Sample_ID, Sample_Name, Fox_Aspect, Month)) %>%
filter(!is.na(Observed)) %>%
group_by(DNA_or_cDNA, Snowpack_Layer) %>%
mutate(Observed = mean(Observed),
Shannon = mean(Shannon),
Simpson = mean(Simpson)) %>%
distinct_all() %>%
kable(digits = 4) %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>%
row_spec(0, bold=TRUE)
## Warning in `==.default`(Month, c("Early July", "Late July")): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning in `!=.default`(Sample_Name, c("30 July_SWS1SE_Slush_R", "30
## July_SWS1SE_Slush_D")): longer object length is not a multiple of shorter object
## length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
|
Observed
|
Shannon
|
Simpson
|
DNA_or_cDNA
|
Snowpack_Layer
|
|
1635.0000
|
4.6591
|
0.9488
|
cDNA
|
MID
|
|
1875.0000
|
5.5902
|
0.9807
|
cDNA
|
TOP
|
|
2250.0000
|
5.2598
|
0.9659
|
DNA
|
TOP
|
|
592.6667
|
3.8378
|
0.9212
|
cDNA
|
SUP ICE
|
|
1586.0000
|
4.6213
|
0.9348
|
DNA
|
MID
|
All measures for Snowpack_Layer in “Early July”, stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed
plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
Month == "Early July" & !(Sample_Name %in% gl_ice_spls)),
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = slayer_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Snowpack Layer") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)
Table with means
estimate_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
Month == "Early July" & !(Sample_Name %in% gl_ice_spls)),
measures = c("Observed","Shannon","Simpson")) %>%
rownames_to_column(var = "Sample_ID") %>%
mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>%
full_join(., mini_metadata, by="Sample_ID") %>%
select(!c(Sample_ID, Sample_Name, Fox_Aspect, Month)) %>%
filter(!is.na(Observed)) %>%
group_by(DNA_or_cDNA, Snowpack_Layer) %>%
mutate(Observed = mean(Observed),
Shannon = mean(Shannon),
Simpson = mean(Simpson)) %>%
distinct_all() %>%
kable(digits = 4) %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>%
row_spec(0, bold=TRUE)
|
Observed
|
Shannon
|
Simpson
|
DNA_or_cDNA
|
Snowpack_Layer
|
|
850.0000
|
3.1169
|
0.8073
|
cDNA
|
MID
|
|
1099.0000
|
3.7314
|
0.8676
|
DNA
|
SUP ICE
|
|
1160.6667
|
3.7859
|
0.8616
|
DNA
|
MID
|
|
592.6667
|
3.8378
|
0.9212
|
cDNA
|
SUP ICE
|
|
1726.5000
|
3.8925
|
0.7601
|
DNA
|
TOP
|
All measures for Snowpack_Layer in Late July, stats included, facet by DNA/cDNA, “GL ICE” removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_R” removed, “Gl ice” samples removed
plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
Month == "Late July" & !(Sample_Name %in% gl_ice_spls)),
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = slayer_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Snowpack Layer") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations

# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)
Table with means
estimate_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
Month == "Late July" & !(Sample_Name %in% gl_ice_spls)),
measures = c("Observed","Shannon","Simpson")) %>%
rownames_to_column(var = "Sample_ID") %>%
mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>%
full_join(., mini_metadata, by="Sample_ID") %>%
select(!c(Sample_ID, Sample_Name, Fox_Aspect, Month)) %>%
filter(!is.na(Observed)) %>%
group_by(DNA_or_cDNA, Snowpack_Layer) %>%
mutate(Observed = mean(Observed),
Shannon = mean(Shannon),
Simpson = mean(Simpson)) %>%
distinct_all() %>%
kable(digits = 4) %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>%
row_spec(0, bold=TRUE)
|
Observed
|
Shannon
|
Simpson
|
DNA_or_cDNA
|
Snowpack_Layer
|
|
952
|
3.2764
|
0.6659
|
cDNA
|
TOP
|
|
2719
|
5.7179
|
0.9837
|
DNA
|
TOP
|
|
453
|
4.0612
|
0.9622
|
cDNA
|
SUP ICE
|
|
1739
|
5.2347
|
0.9838
|
DNA
|
SUP ICE
|
8 - Alpha diversity with rarefaction to the smallest number of reads in any sample
ps.final.rare = rarefy_even_depth(ps.final,
sample.size = min(sample_sums(ps.final)),
rngseed = TRUE, replace = TRUE,
trimOTUs = TRUE, verbose = TRUE)
## `set.seed(TRUE)` was used to initialize repeatable random subsampling.
## Please record this for your records so others can reproduce.
## Try `set.seed(TRUE); .Random.seed` for the full vector
## ...
## 25951OTUs were removed because they are no longer
## present in any sample after random subsampling
## ...
plot_richness(ps.final.rare,
x="Sample_ID", color="Sample_ID",
measures=c("Observed","Shannon","Simpson")) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Sample_ID") +
guides(color=guide_legend(ncol=2)) +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))

plot_richness(ps.final.rare,
x="Month", color="Month",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(shape=16, position=position_jitter(0.2)) +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Month") +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))

plot_richness(ps.final.rare,
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(shape=16, position=position_jitter(0.2)) +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = site_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Fox Aspect") +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))

#now without stats, only shapes
plot_richness(ps.final.rare,
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(aes(shape=Month), size=5, alpha=0.7,
position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Fox Aspect") +
scale_shape_manual(values=nmds_shapes) +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))

plot_richness(ps.final.rare,
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(position=position_jitter(0.2)) +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = slayer_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Snowpack Layer") +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))

#now without stats, only shapes
plot_richness(ps.final.rare,
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(aes(shape=Month), size=5, alpha=0.7,
position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Snowpack Layer") +
scale_shape_manual(values=nmds_shapes) +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))

plot_richness(ps.final.rare,
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(aes(shape=Fox_Aspect), size=5, alpha=0.7,
position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Snowpack Layer") +
scale_shape_manual(values=nmds_shapes) +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
